Importing the Data

# @title
import pandas as pd
from sklearn.model_selection import cross_val_score, TimeSeriesSplit

data = pd.read_csv(
    'https://raw.githubusercontent.com/Garve/datasets/4576d323bf2b66c906d5130d686245ad205505cf/mmm.csv',
    parse_dates=['Date'],
    index_col='Date'
)

X = data.drop(columns=['Sales'])
y = data['Sales']

Parameter Values

# @title
N = 200 # number of observations

# previous marketing mix modeling has given us these parameters
tv_coef = 10000       # α
tv_lags = 4           # ℓ
tv_carryover = 0.5    # λ
tv_saturation = 0.002 # β

radio_coef = 8000
radio_lags = 2
radio_carryover = 0.2
radio_saturation = 0.0001

banners_coef = 14000
banners_lags = 0
banners_carryover = 0.2
banners_saturation = 0.001

Carryover Matrices

# @title
from scipy.linalg import toeplitz
import numpy as np

# Much cleaner!
tv_carryover_matrix = toeplitz(
    tv_carryover**np.arange(N),  # first column
    np.zeros(N)                   # first row (all zeros = lower triangular)
)

radio_carryover_matrix = toeplitz(
    radio_carryover**np.arange(N),
    np.zeros(N)
)

banners_carryover_matrix = np.eye(N)

Prediction

# @title
# Calculate predicted sales per period (corrected formula)
sales_per_period = (
    tv_coef * (1 - np.exp(-tv_saturation * tv_carryover_matrix @ data["TV"]))
    + radio_coef * (1 - np.exp(-radio_saturation * radio_carryover_matrix @ data["Radio"]))
    + banners_coef * (1 - np.exp(-banners_saturation * banners_carryover_matrix @ data["Banners"]))
)

# Total sales comparison
total_actual_sales = y.sum()
total_predicted_sales = sales_per_period.sum()

print(f"Total Actual Sales:    ${total_actual_sales:,.2f}")
print(f"Total Predicted Sales: ${total_predicted_sales:,.2f}")
print(f"Difference:            ${total_predicted_sales - total_actual_sales:,.2f}")
print(f"Percentage Difference: {((total_predicted_sales - total_actual_sales) / total_actual_sales * 100):.2f}%")

Total Actual Sales:    $2,133,628.30
Total Predicted Sales: $3,847,201.38
Difference:            $1,713,573.08
Percentage Difference: 80.31%

Prediction using Optimized Parameters

# @title
from scipy.optimize import minimize
from scipy.linalg import toeplitz
import numpy as np

# Prepare data
X_tv = data["TV"].values
X_radio = data["Radio"].values
X_banners = data["Banners"].values
y_actual = y.values
N = len(y)

def mmm_prediction(params, X_tv, X_radio, X_banners, N):
    """Generate sales prediction given parameters"""
    tv_coef, tv_carryover, tv_saturation = params[0:3]
    radio_coef, radio_carryover, radio_saturation = params[3:6]
    banners_coef, banners_saturation = params[6:8]

    # Create carryover matrices with fixed lags
    tv_lags, radio_lags = 4, 2

    tv_matrix = toeplitz(tv_carryover**np.arange(N), np.zeros(N))
    radio_matrix = toeplitz(radio_carryover**np.arange(N), np.zeros(N))
    banners_matrix = np.eye(N)

    # Calculate contributions
    tv_contribution = tv_coef * (1 - np.exp(-tv_saturation * tv_matrix @ X_tv))
    radio_contribution = radio_coef * (1 - np.exp(-radio_saturation * radio_matrix @ X_radio))
    banners_contribution = banners_coef * (1 - np.exp(-banners_saturation * banners_matrix @ X_banners))

    return tv_contribution + radio_contribution + banners_contribution

def objective(params, X_tv, X_radio, X_banners, y_actual, N):
    """Mean Absolute Percentage Error"""
    y_pred = mmm_prediction(params, X_tv, X_radio, X_banners, N)
    mape = np.mean(np.abs((y_actual - y_pred) / y_actual)) * 100
    return mape

# Initial guess (your original parameters)
initial_params = [
    10000, 0.5, 0.002,      # TV: coef, carryover, saturation
    8000, 0.2, 0.0001,       # Radio: coef, carryover, saturation
    14000, 0.001             # Banners: coef, saturation
]

# Bounds to keep parameters reasonable
bounds = [
    (0, 50000),      # tv_coef
    (0, 0.99),       # tv_carryover
    (0.0001, 0.01),  # tv_saturation
    (0, 50000),      # radio_coef
    (0, 0.99),       # radio_carryover
    (0.0001, 0.01),  # radio_saturation
    (0, 50000),      # banners_coef
    (0.0001, 0.01),  # banners_saturation
]

print("Optimizing parameters...")
result = minimize(
    objective,
    initial_params,
    args=(X_tv, X_radio, X_banners, y_actual, N),
    method='L-BFGS-B',
    bounds=bounds,
    options={'maxiter': 1000}
)

# Extract optimized parameters
opt_params = result.x
print(f"\nOptimization successful: {result.success}")
print(f"Final MAPE: {result.fun:.2f}%")
print("\nOptimized Parameters:")
print(f"TV:      coef={opt_params[0]:.2f}, carryover={opt_params[1]:.4f}, saturation={opt_params[2]:.6f}")
print(f"Radio:   coef={opt_params[3]:.2f}, carryover={opt_params[4]:.4f}, saturation={opt_params[5]:.6f}")
print(f"Banners: coef={opt_params[6]:.2f}, saturation={opt_params[7]:.6f}")

# Compare results
y_pred_optimized = mmm_prediction(opt_params, X_tv, X_radio, X_banners, N)

print("\n=== Total Sales Comparison ===")
print(f"Actual:              ${y_actual.sum():,.2f}")
print(f"Original Model:      ${sales_per_period.sum():,.2f}")
print(f"Optimized Model:     ${y_pred_optimized.sum():,.2f}")

Optimizing parameters...

Optimization successful: True
Final MAPE: 12.19%

Optimized Parameters:
TV:      coef=10000.00, carryover=0.7257, saturation=0.000100
Radio:   coef=8000.00, carryover=0.2909, saturation=0.000104
Banners: coef=14000.00, saturation=0.000145

=== Total Sales Comparison ===
Actual:              $2,133,628.30
Original Model:      $3,847,201.38
Optimized Model:     $2,085,086.69

# @title
import matplotlib.pyplot as plt

fig, axes = plt.subplots(2, 1, figsize=(14, 8))

# Time series comparison
axes[0].plot(data.index, y_actual, label='Actual', linewidth=2)
axes[0].plot(data.index, y_pred_optimized, label='Predicted', linewidth=2, alpha=0.7)
axes[0].set_title('Actual vs Predicted Sales Over Time')
axes[0].legend()
axes[0].set_ylabel('Sales')

# Residuals
residuals = y_actual - y_pred_optimized
axes[1].plot(data.index, residuals)
axes[1].axhline(y=0, color='r', linestyle='--')
axes[1].set_title('Residuals')
axes[1].set_ylabel('Residual')
plt.tight_layout()
plt.show()

Contributions by channel

# @title
# Recreate the carryover matrices with optimized parameters
tv_matrix = toeplitz(opt_params[1]**np.arange(N), np.zeros(N))
radio_matrix = toeplitz(opt_params[4]**np.arange(N), np.zeros(N))
banners_matrix = np.eye(N)

# Calculate individual channel contributions
tv_contribution = opt_params[0] * (1 - np.exp(-opt_params[2] * tv_matrix @ X_tv))
radio_contribution = opt_params[3] * (1 - np.exp(-opt_params[5] * radio_matrix @ X_radio))
banners_contribution = opt_params[6] * (1 - np.exp(-opt_params[7] * banners_matrix @ X_banners))

print("Total Sales Contribution by Channel:")
print(f"TV:      ${tv_contribution.sum():,.2f} ({tv_contribution.sum()/y_pred_optimized.sum()*100:.1f}%)")
print(f"Radio:   ${radio_contribution.sum():,.2f} ({radio_contribution.sum()/y_pred_optimized.sum()*100:.1f}%)")
print(f"Banners: ${banners_contribution.sum():,.2f} ({banners_contribution.sum()/y_pred_optimized.sum()*100:.1f}%)")

Total Sales Contribution by Channel:
TV:      $1,146,349.83 (55.0%)
Radio:   $403,195.47 (19.3%)
Banners: $535,541.40 (25.7%)

Goodness of Fit

# @title
# Calculate R-squared
ss_res = np.sum((y_actual - y_pred_optimized)**2)  # Residual sum of squares
ss_tot = np.sum((y_actual - y_actual.mean())**2)   # Total sum of squares
r_squared = 1 - (ss_res / ss_tot)

print("\n=== Model Performance Metrics ===")
print(f"R-squared:           {r_squared:.4f}")
print(f"Adjusted R-squared:  {1 - (1-r_squared)*(N-1)/(N-8-1):.4f}")  # 8 parameters
print(f"MAPE:                {result.fun:.2f}%")
print(f"RMSE:                ${np.sqrt(ss_res/N):,.2f}")

# Alternative: using sklearn
from sklearn.metrics import r2_score, mean_absolute_percentage_error, mean_squared_error

print("\n=== Using sklearn ===")
print(f"R²:                  {r2_score(y_actual, y_pred_optimized):.4f}")
print(f"MAPE:                {mean_absolute_percentage_error(y_actual, y_pred_optimized)*100:.2f}%")
print(f"RMSE:                ${np.sqrt(mean_squared_error(y_actual, y_pred_optimized)):,.2f}")

=== Model Performance Metrics ===
R-squared:           0.7101
Adjusted R-squared:  0.6979
MAPE:                12.19%
RMSE:                $1,450.54

=== Using sklearn ===
R²:                  0.7101
MAPE:                12.19%
RMSE:                $1,450.54

Optimal Budget Allocation Using CVXPY

# @title
import cvxpy as cp
import numpy as np
from scipy.linalg import toeplitz

# Extract optimized parameters
tv_coef = opt_params[0]
tv_carryover = opt_params[1]
tv_saturation = opt_params[2]

radio_coef = opt_params[3]
radio_carryover = opt_params[4]
radio_saturation = opt_params[5]

banners_coef = opt_params[6]
banners_saturation = opt_params[7]

# Create carryover matrices with optimized parameters
tv_carryover_matrix = toeplitz(tv_carryover**np.arange(N), np.zeros(N))
radio_carryover_matrix = toeplitz(radio_carryover**np.arange(N), np.zeros(N))
banners_carryover_matrix = np.eye(N)

# Get original total budget
original_total_spends = data[["TV", "Radio", "Banners"]].sum().sum()

# Declare CVXPY variables (N=200 per channel)
tv = cp.Variable(N)
radio = cp.Variable(N)
banners = cp.Variable(N)

# Constraints
constraints = [
    tv >= 0,
    radio >= 0,
    banners >= 0,
    cp.sum(tv + radio + banners) <= original_total_spends,
]

# Objective function using optimized parameters
objective = cp.Maximize(
    tv_coef * cp.sum(1 - cp.exp(-tv_saturation * tv_carryover_matrix @ tv))
    + radio_coef * cp.sum(1 - cp.exp(-radio_saturation * radio_carryover_matrix @ radio))
    + banners_coef * cp.sum(1 - cp.exp(-banners_saturation * banners_carryover_matrix @ banners))
)

# Create and solve problem (CVXPY auto-selects solver)
problem = cp.Problem(objective, constraints)
problem.solve(verbose=True)

# Display results
print("\n" + "="*60)
print("OPTIMIZATION RESULTS")
print("="*60)
print(f"Solver used: {problem.solver_stats.solver_name}")
print(f"Status: {problem.status}")
print(f"\nOptimal Total Sales: ${problem.value:,.2f}")
print(f"\n--- Budget Allocation ---")
print(f"Original Total Budget:  ${original_total_spends:,.2f}")
print(f"Optimized Total Budget: ${(tv.value.sum() + radio.value.sum() + banners.value.sum()):,.2f}")
print(f"\n--- Spend per Channel ---")
print(f"TV:      ${tv.value.sum():,.2f} (Original: ${data['TV'].sum():,.2f})")
print(f"Radio:   ${radio.value.sum():,.2f} (Original: ${data['Radio'].sum():,.2f})")
print(f"Banners: ${banners.value.sum():,.2f} (Original: ${data['Banners'].sum():,.2f})")
print(f"\n--- Sales Improvement ---")
print(f"Actual Sales (from data):    ${y_actual.sum():,.2f}")
print(f"Optimized Sales:             ${problem.value:,.2f}")
print(f"Improvement:                 ${problem.value - y_actual.sum():,.2f}")
print(f"Improvement %:               {((problem.value - y_actual.sum()) / y_actual.sum() * 100):.2f}%")

(CVXPY) Jan 17 01:28:50 PM: Your problem has 600 variables, 601 constraints, and 0 parameters.
(CVXPY) Jan 17 01:28:50 PM: It is compliant with the following grammars: DCP, DQCP
(CVXPY) Jan 17 01:28:50 PM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)
(CVXPY) Jan 17 01:28:50 PM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.
(CVXPY) Jan 17 01:28:50 PM: Your problem is compiled with the CPP canonicalization backend.
(CVXPY) Jan 17 01:28:50 PM: Compiling problem (target solver=CLARABEL).
(CVXPY) Jan 17 01:28:50 PM: Reduction chain: FlipObjective -> Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> CLARABEL
(CVXPY) Jan 17 01:28:50 PM: Applying reduction FlipObjective
(CVXPY) Jan 17 01:28:50 PM: Applying reduction Dcp2Cone
(CVXPY) Jan 17 01:28:50 PM: Applying reduction CvxAttr2Constr
(CVXPY) Jan 17 01:28:50 PM: Applying reduction ConeMatrixStuffing
(CVXPY) Jan 17 01:28:50 PM: Applying reduction CLARABEL
(CVXPY) Jan 17 01:28:50 PM: Finished problem compilation (took 1.177e-01 seconds).
(CVXPY) Jan 17 01:28:50 PM: Invoking solver CLARABEL  to obtain a solution.

===============================================================================
                                     CVXPY                                     
                                     v1.6.7                                    
===============================================================================
-------------------------------------------------------------------------------
                                  Compilation                                  
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
                                Numerical solver                               
-------------------------------------------------------------------------------
-------------------------------------------------------------
           Clarabel.rs v0.11.1  -  Clever Acronym                

                   (c) Paul Goulart                          
                University of Oxford, 2022                   
-------------------------------------------------------------

problem:
  variables     = 1200
  constraints   = 2401
  nnz(P)        = 0
  nnz(A)        = 42200
  cones (total) = 601
    : Nonnegative = 1,  numel = 601
    : Exponential = 600,  numel = (3,3,3,3,...,3)

settings:
  linear algebra: direct / faer, precision: 64 bit (2 threads)
  max iter = 200, time limit = Inf,  max step = 0.990
  tol_feas = 1.0e-8, tol_gap_abs = 1.0e-8, tol_gap_rel = 1.0e-8,
  static reg : on, ϵ1 = 1.0e-8, ϵ2 = 4.9e-32
  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-7
  iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
               max iter = 10, stop ratio = 5.0
  equilibrate: on, min_scale = 1.0e-4, max_scale = 1.0e4
               max iter = 10

iter    pcost        dcost       gap       pres      dres      k/t        μ       step      
---------------------------------------------------------------------------------------------
  0  +0.0000e+00  -2.3105e+06  2.31e+06  1.00e+00  1.34e+00  1.00e+00  1.00e+00   ------   
  1  +1.2308e+04  -2.2963e+06  1.88e+02  9.99e-01  1.04e+00  4.08e+00  2.45e-01  7.92e-01  
  2  +3.0658e+05  -1.9551e+06  7.38e+00  9.77e-01  1.26e+00  7.09e+01  1.38e-02  9.90e-01  
  3  +1.1475e+06  -9.5952e+05  2.20e+00  9.06e-01  1.09e+00  2.67e+02  3.29e-03  7.71e-01  
  4  +2.4606e+06  +7.6837e+05  2.20e+00  7.17e-01  6.44e-01  5.87e+02  1.03e-03  7.12e-01  
  5  +3.4285e+06  +2.6422e+06  2.98e-01  3.23e-01  1.83e-01  6.33e+02  2.44e-04  7.99e-01  
  6  +3.7013e+06  +3.3954e+06  9.01e-02  1.23e-01  5.85e-02  2.51e+02  7.63e-05  7.92e-01  
  7  +3.8622e+06  +3.7795e+06  2.19e-02  3.28e-02  1.46e-02  6.98e+01  1.88e-05  7.73e-01  
  8  +3.9037e+06  +3.8820e+06  5.58e-03  8.55e-03  3.72e-03  1.83e+01  4.80e-06  7.64e-01  
  9  +3.9143e+06  +3.9094e+06  1.27e-03  1.96e-03  8.48e-04  4.22e+00  1.09e-06  7.92e-01  
 10  +3.9165e+06  +3.9154e+06  2.82e-04  4.34e-04  1.88e-04  9.36e-01  2.42e-07  7.92e-01  
 11  +3.9170e+06  +3.9167e+06  6.12e-05  9.43e-05  4.07e-05  2.04e-01  5.25e-08  7.92e-01  
 12  +3.9170e+06  +3.9170e+06  1.32e-05  2.03e-05  8.78e-06  4.39e-02  1.13e-08  7.92e-01  
 13  +3.9171e+06  +3.9170e+06  2.85e-06  4.38e-06  1.89e-06  9.47e-03  2.44e-09  7.91e-01  
 14  +3.9171e+06  +3.9171e+06  6.30e-07  9.71e-07  4.20e-07  2.10e-03  5.41e-10  7.85e-01  
 15  +3.9171e+06  +3.9171e+06  3.78e-07  5.82e-07  2.52e-07  1.26e-03  3.24e-10  5.07e-01  
 16  +3.9171e+06  +3.9171e+06  3.51e-07  5.41e-07  2.34e-07  1.17e-03  3.01e-10  2.08e-01  

(CVXPY) Jan 17 01:28:51 PM: Problem status: optimal
(CVXPY) Jan 17 01:28:51 PM: Optimal value: 2.483e+06
(CVXPY) Jan 17 01:28:51 PM: Compilation took 1.177e-01 seconds
(CVXPY) Jan 17 01:28:51 PM: Solver (including time spent in interface) took 6.037e-01 seconds

 17  +3.9171e+06  +3.9171e+06  3.38e-07  5.21e-07  2.25e-07  1.13e-03  2.90e-10  1.06e-01  
 18  +3.9171e+06  +3.9171e+06  3.38e-07  5.21e-07  2.25e-07  1.13e-03  2.90e-10  0.00e+00  
 19  +3.9171e+06  +3.9171e+06  1.33e-07  2.06e-07  8.88e-08  4.42e-04  1.16e-10  6.34e-01  
 20  +3.9171e+06  +3.9171e+06  2.86e-08  4.41e-08  1.90e-08  9.46e-05  2.49e-11  7.92e-01  
 21  +3.9171e+06  +3.9171e+06  1.13e-08  1.75e-08  7.55e-09  3.76e-05  9.85e-12  6.34e-01  
 22  +3.9171e+06  +3.9171e+06  2.43e-09  3.74e-09  1.62e-09  8.06e-06  2.12e-12  7.92e-01  
---------------------------------------------------------------------------------------------
Terminated with status = Solved
solve time = 582.139145ms
-------------------------------------------------------------------------------
                                    Summary                                    
-------------------------------------------------------------------------------

============================================================
OPTIMIZATION RESULTS
============================================================
Solver used: CLARABEL
Status: optimal

Optimal Total Sales: $2,482,938.10

--- Budget Allocation ---
Original Total Budget:  $1,336,103.05
Optimized Total Budget: $1,336,103.04

--- Spend per Channel ---
TV:      $609,141.90 (Original: $589,241.53)
Radio:   $0.08 (Original: $442,717.01)
Banners: $726,961.06 (Original: $304,144.51)

--- Sales Improvement ---
Actual Sales (from data):    $2,133,628.30
Optimized Sales:             $2,482,938.10
Improvement:                 $349,309.80
Improvement %:               16.37%

Budget Allocation Optimization with Scipy Minimize (SLSQP)

# @title
from scipy.optimize import minimize
from functools import partial
import numpy as np
from scipy.linalg import toeplitz

# Extract optimized parameters
tv_coef = opt_params[0]
tv_carryover = opt_params[1]
tv_saturation = opt_params[2]

radio_coef = opt_params[3]
radio_carryover = opt_params[4]
radio_saturation = opt_params[5]

banners_coef = opt_params[6]
banners_saturation = opt_params[7]

# Create carryover matrices
N = 200
tv_carryover_matrix = toeplitz(tv_carryover**np.arange(N), np.zeros(N))
radio_carryover_matrix = toeplitz(radio_carryover**np.arange(N), np.zeros(N))
banners_carryover_matrix = np.eye(N)

# Get original budget and channel totals
original_total_budget = data[["TV", "Radio", "Banners"]].sum().sum()
original_channel_spends = data[["TV", "Radio", "Banners"]].sum().values

# Initial guess - original total spends per channel
initial_channel_totals = original_channel_spends

# Bounds - all channel spends must be >= 0
bounds = [(0, None), (0, None), (0, None)]

# Objective function to minimize (we'll negate sales since we want to maximize)
def mmm_objective_by_channel(channel_totals, tv_coef, tv_sat, tv_matrix,
                              radio_coef, radio_sat, radio_matrix,
                              banners_coef, banners_sat, banners_matrix,
                              original_spends, N):
    """
    Calculate negative total sales for optimization.
    channel_totals: [tv_total, radio_total, banners_total]
    Distributes each channel's total proportionally to original daily spends
    """
    tv_total, radio_total, banners_total = channel_totals

    # Distribute budget proportionally to original daily pattern
    tv_spend = original_spends["TV"].values * (tv_total / original_spends["TV"].sum())
    radio_spend = original_spends["Radio"].values * (radio_total / original_spends["Radio"].sum())
    banners_spend = original_spends["Banners"].values * (banners_total / original_spends["Banners"].sum())

    # Calculate contributions with adstock and saturation
    tv_contribution = tv_coef * np.sum(1 - np.exp(-tv_sat * tv_matrix @ tv_spend))
    radio_contribution = radio_coef * np.sum(1 - np.exp(-radio_sat * radio_matrix @ radio_spend))
    banners_contribution = banners_coef * np.sum(1 - np.exp(-banners_sat * banners_matrix @ banners_spend))

    total_sales = tv_contribution + radio_contribution + banners_contribution

    # Return negative because we're minimizing
    return -total_sales

# Budget constraint (equality)
def budget_constraint(channel_totals, budget):
    """Equality constraint: total spend must equal budget"""
    return np.sum(channel_totals) - budget

# Create partial function with fixed parameters
partial_objective = partial(
    mmm_objective_by_channel,
    tv_coef=tv_coef,
    tv_sat=tv_saturation,
    tv_matrix=tv_carryover_matrix,
    radio_coef=radio_coef,
    radio_sat=radio_saturation,
    radio_matrix=radio_carryover_matrix,
    banners_coef=banners_coef,
    banners_sat=banners_saturation,
    banners_matrix=banners_carryover_matrix,
    original_spends=data,
    N=N
)

# Solve optimization
print("Starting optimization with scipy.optimize.minimize (SLSQP)...")
solution = minimize(
    fun=partial_objective,
    x0=initial_channel_totals,
    bounds=bounds,
    method="SLSQP",
    jac="3-point",
    options={
        "maxiter": 1000,
        "disp": True,
        "ftol": 1e-9
    },
    constraints={
        "type": "eq",
        "fun": budget_constraint,
        "args": (original_total_budget,)
    }
)

# Extract optimized channel totals
tv_total_optimized, radio_total_optimized, banners_total_optimized = solution.x

# Display results
print("\n" + "="*60)
print("SCIPY OPTIMIZATION RESULTS (SLSQP)")
print("="*60)
print(f"Status: {solution.message}")
print(f"Success: {solution.success}")
print(f"\nOptimal Total Sales: ${-solution.fun:,.2f}")  # Negative because we minimized
print(f"\n--- Budget Allocation ---")
print(f"Original Total Budget:  ${original_total_budget:,.2f}")
print(f"Optimized Total Budget: ${solution.x.sum():,.2f}")
print(f"\n--- Spend per Channel ---")
print(f"TV:      ${tv_total_optimized:,.2f} (Original: ${data['TV'].sum():,.2f})")
print(f"Radio:   ${radio_total_optimized:,.2f} (Original: ${data['Radio'].sum():,.2f})")
print(f"Banners: ${banners_total_optimized:,.2f} (Original: ${data['Banners'].sum():,.2f})")
print(f"\n--- Sales Improvement ---")
print(f"Actual Sales (from data):    ${y_actual.sum():,.2f}")
print(f"Optimized Sales (SLSQP):     ${-solution.fun:,.2f}")
print(f"Improvement:                 ${-solution.fun - y_actual.sum():,.2f}")
print(f"Improvement %:               {((-solution.fun - y_actual.sum()) / y_actual.sum() * 100):.2f}%")

# Compare with CVXPY results
print(f"\n--- Comparison with CVXPY ---")
print(f"CVXPY Sales:  ${problem.value:,.2f}")
print(f"SLSQP Sales:  ${-solution.fun:,.2f}")
print(f"Difference:   ${abs(problem.value - (-solution.fun)):,.2f}")

Starting optimization with scipy.optimize.minimize (SLSQP)...
Optimization terminated successfully    (Exit mode 0)
            Current function value: -2237678.352585092
            Iterations: 32
            Function evaluations: 224
            Gradient evaluations: 32

============================================================
SCIPY OPTIMIZATION RESULTS (SLSQP)
============================================================
Status: Optimization terminated successfully
Success: True

Optimal Total Sales: $2,237,678.35

--- Budget Allocation ---
Original Total Budget:  $1,336,103.05
Optimized Total Budget: $1,336,103.05

--- Spend per Channel ---
TV:      $551,512.71 (Original: $589,241.53)
Radio:   $93,956.04 (Original: $442,717.01)
Banners: $690,634.30 (Original: $304,144.51)

--- Sales Improvement ---
Actual Sales (from data):    $2,133,628.30
Optimized Sales (SLSQP):     $2,237,678.35
Improvement:                 $104,050.05
Improvement %:               4.88%

--- Comparison with CVXPY ---
CVXPY Sales:  $2,482,938.10
SLSQP Sales:  $2,237,678.35
Difference:   $245,259.74

Comparison of Budget Allocation Methods

# @title
import matplotlib.pyplot as plt
import numpy as np

# Prepare data for comparison
methods = ['Original', 'CVXPY', 'SLSQP']

# Budget allocations
tv_spends = [data['TV'].sum(), tv.value.sum(), tv_total_optimized]
radio_spends = [data['Radio'].sum(), radio.value.sum(), radio_total_optimized]
banners_spends = [data['Banners'].sum(), banners.value.sum(), banners_total_optimized]

# Total sales
sales_values = [y_actual.sum(), problem.value, -solution.fun]

# Create comparison plots
fig, axes = plt.subplots(2, 2, figsize=(15, 10), dpi=150)

# 1. Budget Allocation Comparison (Stacked Bar)
x = np.arange(len(methods))
width = 0.5

axes[0, 0].bar(x, tv_spends, width, label='TV', color='#1f77b4')
axes[0, 0].bar(x, radio_spends, width, bottom=tv_spends, label='Radio', color='#ff7f0e')
axes[0, 0].bar(x, banners_spends, width, bottom=np.array(tv_spends)+np.array(radio_spends),
               label='Banners', color='#2ca02c')
axes[0, 0].set_ylabel('Budget ($)', fontsize=11)
axes[0, 0].set_title('Budget Allocation by Method', fontsize=12, fontweight='bold')
axes[0, 0].set_xticks(x)
axes[0, 0].set_xticklabels(methods)
axes[0, 0].legend()
axes[0, 0].grid(axis='y', alpha=0.3)

# 2. Total Sales Comparison (Bar)
colors = ['#666666', '#d62728', '#9467bd']
bars = axes[0, 1].bar(methods, sales_values, color=colors, alpha=0.8, edgecolor='black')
axes[0, 1].set_ylabel('Total Sales ($)', fontsize=11)
axes[0, 1].set_title('Total Sales by Method', fontsize=12, fontweight='bold')
axes[0, 1].grid(axis='y', alpha=0.3)

# Add value labels on bars
for bar, val in zip(bars, sales_values):
    height = bar.get_height()
    axes[0, 1].text(bar.get_x() + bar.get_width()/2., height,
                    f'${val:,.0f}',
                    ha='center', va='bottom', fontsize=9)

# 3. Channel Budget Comparison (Grouped Bar)
x = np.arange(len(methods))
width = 0.25

axes[1, 0].bar(x - width, tv_spends, width, label='TV', color='#1f77b4', alpha=0.8)
axes[1, 0].bar(x, radio_spends, width, label='Radio', color='#ff7f0e', alpha=0.8)
axes[1, 0].bar(x + width, banners_spends, width, label='Banners', color='#2ca02c', alpha=0.8)

axes[1, 0].set_ylabel('Spend ($)', fontsize=11)
axes[1, 0].set_title('Channel Spend Comparison', fontsize=12, fontweight='bold')
axes[1, 0].set_xticks(x)
axes[1, 0].set_xticklabels(methods)
axes[1, 0].legend()
axes[1, 0].grid(axis='y', alpha=0.3)

# 4. Sales Improvement Percentage (Bar)
improvement_cvxpy = ((problem.value - y_actual.sum()) / y_actual.sum()) * 100
improvement_slsqp = ((-solution.fun - y_actual.sum()) / y_actual.sum()) * 100

improvements = [0, improvement_cvxpy, improvement_slsqp]
colors_imp = ['#666666', '#d62728', '#9467bd']
bars = axes[1, 1].bar(methods, improvements, color=colors_imp, alpha=0.8, edgecolor='black')
axes[1, 1].set_ylabel('Sales Improvement (%)', fontsize=11)
axes[1, 1].set_title('Sales Improvement vs Original', fontsize=12, fontweight='bold')
axes[1, 1].axhline(y=0, color='black', linestyle='-', linewidth=0.8)
axes[1, 1].grid(axis='y', alpha=0.3)

# Add value labels on bars
for bar, val in zip(bars, improvements):
    height = bar.get_height()
    axes[1, 1].text(bar.get_x() + bar.get_width()/2., height,
                    f'{val:.2f}%',
                    ha='center', va='bottom' if val > 0 else 'top', fontsize=9)

plt.tight_layout()
plt.show()

# Print detailed comparison table
print("\n" + "="*80)
print("DETAILED COMPARISON: CVXPY vs SLSQP")
print("="*80)
print(f"{'Metric':<30} {'CVXPY':<20} {'SLSQP':<20} {'Difference':<15}")
print("-"*80)
print(f"{'Total Sales':<30} ${problem.value:>15,.2f} ${-solution.fun:>15,.2f} ${abs(problem.value - (-solution.fun)):>12,.2f}")
print(f"{'TV Spend':<30} ${tv.value.sum():>15,.2f} ${tv_total_optimized:>15,.2f} ${abs(tv.value.sum() - tv_total_optimized):>12,.2f}")
print(f"{'Radio Spend':<30} ${radio.value.sum():>15,.2f} ${radio_total_optimized:>15,.2f} ${abs(radio.value.sum() - radio_total_optimized):>12,.2f}")
print(f"{'Banners Spend':<30} ${banners.value.sum():>15,.2f} ${banners_total_optimized:>15,.2f} ${abs(banners.value.sum() - banners_total_optimized):>12,.2f}")
print(f"{'Sales Improvement %':<30} {improvement_cvxpy:>15.2f}% {improvement_slsqp:>15.2f}% {abs(improvement_cvxpy - improvement_slsqp):>12.2f}%")
print("="*80)

================================================================================
DETAILED COMPARISON: CVXPY vs SLSQP
================================================================================
Metric                         CVXPY                SLSQP                Difference     
--------------------------------------------------------------------------------
Total Sales                    $   2,482,938.10 $   2,237,678.35 $  245,259.74
TV Spend                       $     609,141.90 $     551,512.71 $   57,629.20
Radio Spend                    $           0.08 $      93,956.04 $   93,955.96
Banners Spend                  $     726,961.06 $     690,634.30 $   36,326.76
Sales Improvement %                      16.37%            4.88%        11.49%
================================================================================

CVXPY Constraints:

Total: 601 constraints Budget constraint type: Inequality (≤) - can spend up to the budget Variable count: 600 variables (200 time periods × 3 channels)

scipy.optimize.minimize Constraints:

Total: 3 bound constraints + 1 equality constraint = 4 constraints Budget constraint type: Equality (=) - must spend exactly the budget Variable count: 3 variables (total per channel)

Budget Optimization Analysis | Marketing Mix Modeling | 2026

Built with Jupyter Notebook | Optimized for GitHub Pages